#ifndef CLASS_ICS
#define CLASS_ICS

#ifdef DOUBLE
#define REAL double
#define INT  lint
#define UINT luint
#define __SFT 52
#endif

#ifdef FLOAT
#define REAL float
#define INT  int
#define UINT uint
#define __SFT 23
#endif


#define SQR(x) ((x)*(x))
#define MOD3(x) ( sqrt((x[0])*(x[0]) + (x[1])*(x[1]) + (x[2])*(x[2])) )

#include <xmmintrin.h>
#include <stdio.h>
#include <dirent.h>
#include <iostream>


#include "Timer.h"
#include "Halo.h"
#include "Disk.h"
#include "Vector.h"
#include "AG.h"
#include "DG.h"
#include "GadgetFormat1IO.h"
#include "FileListFinder.h"
#include "FileReader.h"



typedef unsigned int uint;
typedef long int lint;
typedef long unsigned int luint;

using namespace std;


class ICs {
	public:
	
		
	double lU, mU, vU; 
	
	REAL c;				
	REAL v200;

	uint halo_n;
	REAL halo_rmax;
	
	uint   disk_n;
	REAL disk_mu;
	REAL disk_r0;
	REAL disk_rmax;
	REAL disk_z0;
	REAL disk_zmax;
		
	uint   bulge_n;
	REAL bulge_mu;
	REAL bulge_r0;
	REAL bulge_rmax;

	
	// integrator
	REAL eps; // softening
	int ns; // steps per period
	int np; // periods per orbit
	
	
	// optimization
 	REAL dvk_mu; // incriment vector magnitude as frac of v_esc  
	REAL dvh_mu; // dv_move max as frac of vj  
	REAL fLz;
	int    nsteps; // number of iterations
	int	 stype; //type of S function
	int	 sft;
	REAL dh_min;
	REAL dh_max;
	int step_out;
	REAL h_out;
	
	char   curfolder[256]; // current folder
	char    ssfolder[256]; // snapshot folder
	

	// track dump
	bool trackF;
	int track_freq; 
	int tracks_freq;
	char trackfolder[256]; // tracks folder
	
	
	
	// FIXED VARs:
	static const double pi = 3.14159265358979323846;
	static const double GRAVITY = 6.672e-8;
	static const double SOLAR_MASS = 1.989e33;
	static const double CM_PER_MPC = 3.085678e24;    
	static const double HUBBLE = 1.0;     // Hubble constant in 100km/sec/Mpc 
	static const int GB = 1024*1024*1024;
	
	// DERIVATIVE VARs
	int n;
	double tU;
	double H0;
	double G;
	
	REAL r200;
	REAL m200;
	
	//REAL M0;
	
	int wn;
	int step;
	int i_out;

	
	int nn; //number of cells intersected by orbit
	
	REAL h, dh;
	
	
	// arrays
	REAL *m0;
	REAL *m;	
	REAL *ms;
	
	
	float *w, *wa;
	float *track;
	
	bool resumeF;
	
	
	// COMPLEX VARs
	
	P *p;
	
	union DoubleLuintUnion {
		REAL f;
		luint i;	
	};
	
	Halo halo;
	Disk disk;
	
	AG ag;
	DG dg;
	
	
	class SG {
		public:
			int n;
			REAL r1, r2;
			SG(): n(8), r1(1), r2(10000) {}
	} sg;
	
	class VG {
		public:
			int n;
			VG(): n(20){}
	} vg;
		
	// MEMBER FUNCTIONS
		
	ICs() :
	p(0),
	wa(0),
	step(0),
	h(0),
	fLz(0),
	i_out(0),
	resumeF(false)
	{
		
		// defaults
		mU = 1e10*1.989e33;	// 10^10 solar masses [g] 
		lU = CM_PER_MPC/1e3;	// 1 kpc  [cm]
		vU = 1e5;	// 1 km/sec [cm]
		
		v200 = 160; // km/s			
		c = 10;
	
		halo_rmax = 100000;
		disk_mu = 0.04;
		bulge_mu = 0;

		
		step_out = 100;
		h_out = 0;
		stype = 2;
		eps = 0.01;
		
		ns = 250;
		np = 75;
		
		nsteps = 100000;
		sft = 2;
		
		dvk_mu = 0.002;
		dvh_mu = 0.002;
	
		dh_min = 1.0/64;
		dh_max = 1.0/2;
		
		
		trackF = false;
		track_freq = 1;
		tracks_freq = 25; 
		
		sprintf(curfolder, "./");
		sprintf(ssfolder, "out");
		sprintf(trackfolder, "tracks");
		
		printf("ICs::ICs()\n");
		
		update_derivatives();
		

	}
	
	void update_derivatives() {
		
		
		tU = lU/vU;
					
		H0 = HUBBLE * 100 * 1e5 / CM_PER_MPC / vU * lU;
		G  = GRAVITY / pow(lU, 3) * mU * pow(tU, 2);
		
		r200 = v200 / (10*H0);
		m200 = pow(v200, 3) / (10*G*H0);
		
		halo.n = halo_n;
		halo.r0 = r200/c * sqrt(2 * (log(1 + c) - c / (1 + c)));
		halo.m = m200*(1-disk_mu-bulge_mu);
		
		
		
		halo.rmax = halo_rmax;
		halo.G = G;

		disk.n = disk_n;
		disk.m = m200*disk_mu;
		disk.r0 = disk_r0;
		disk.z0 = disk_z0;
		disk.rmax = disk_rmax;
		disk.zmax = disk_zmax;
		disk.G = G;
		
		n = halo_n;
		
		
		ag.update_derivatives();
		dg.update_derivatives();
		
	}

		#include ICS_CALC_ORB

		#include ICS_FIND_VEL 
		
		// save array as a float
		void save_arr(string filename, REAL *arr, int size){
			
			printf("save_arr('%s', 0x%x, %i)...", filename.c_str(), *(INT*)arr, size); fflush(stdout);
			
			float *arr_float = new float[size];
			for (int i=0; i<size; i++) arr_float[i] = arr[i];

			FILE *fp = _fopen(filename, "wb");
			fwrite(arr_float, sizeof(float), size, fp);	
			fclose(fp);
			
			delete [] arr_float;
		
			printf("done\n");
		}
		
		
		/*
		void save_dm(string filename){
			
			printf("save_drho_float('%s')...", filename.c_str()); fflush(stdout);
			int dm_size = dg.nr*dg.nz;
			
			float *dm = new float[dm_size];
			for (int i=0; i<dm_size; i++) dm[i] = ms[i]-m0[i];

			FILE *fp = _fopen(filename, "wb");
			fwrite(drho_float, sizeof(float), dm_size, fp);	
			fclose(fp);
			
			delete [] dm;
		
			printf("done\n");
		}
		*/		
		void calc_w(float *w) {
			
			printf("calc_w()..."); fflush(stdout);
		
			wn = sg.n*vg.n;
			
			REAL snrlogr2r1 = sg.n/log(sg.r2/sg.r1); 
			REAL slogr1 = log(sg.r1); 
			
			memset(w, 0, wn*sizeof(*w) );
			
			for (int i=0; i<n; i++) {
			
				REAL r = sqrt(	SQR(p[i].Pos[0]) + 
										SQR(p[i].Pos[1]) + 
										SQR(p[i].Pos[2]) );
				
				int ir =  snrlogr2r1*(log(r)-slogr1);
			
				if (ir < 0  || sg.n <= ir) continue;
				
				REAL v = (p[i].Vel[0]*p[i].Pos[0] + 
								p[i].Vel[1]*p[i].Pos[1] + 
								p[i].Vel[2]*p[i].Pos[2])/r;
				
				REAL rc = exp(slogr1 + (ir+0.5)*log(sg.r2/sg.r1)/sg.n );
				
				REAL vr_max  = halo.vesc(rc);
				
				REAL dv = 2*vr_max/vg.n;
				
				int vid = int(floor(v/dv) ) + vg.n/2;
				
				if (vid < 0 || vg.n <= vid) continue;
				
				int id = ir*vg.n + vid;
				
				w[id]++;
				
			}

			for (int i=0; i<sg.n; i++){
				REAL ws(0); 
				for (int k=0; k<vg.n; k++) ws+= w[i*vg.n + k];
				for (int k=0; k<vg.n; k++) w[i*vg.n + k]/=ws ;
			}
			
			printf("done\n");
		}
		
		
		void save_w(float *w, string filename){
			
			printf("save_w_float('%s')...", filename.c_str()); fflush(stdout);
			
			FILE *fp = _fopen(filename, "wb");
			fwrite(w, sizeof(*w), sg.n*vg.n, fp);	
			fclose(fp);
			
			printf("done\n");
		}

		/*
		void save_ms(string filename){
			
			printf("save_rhos_float('%s')...", filename.c_str()); fflush(stdout);
			
			int ms_size = dg.nr*dg.nz;
			
			float *ms_float = new float[ms_size];
			for (int i=0; i<rhos_size; i++) ms_float[i] = ms[i];

			FILE *fp = _fopen(filename, "wb");
			fwrite(rhos_float, sizeof(float), rhos_size, fp);	
			fclose(fp);

			delete [] rhos_float;
			
			printf("done\n");
		}*/
		
		void save_vel(string filename){
				
			printf("save_vel_float('%s')...", filename.c_str()); fflush(stdout);
			
			FILE *fp = _fopen(filename, "wb");
			
			for (int i=0; i<n; i++) {
				float vel[] = { p[i].Vel[0], p[i].Vel[1], p[i].Vel[2], 1 };
				fwrite(vel, sizeof(*w), 4, fp);	
			}
			fclose(fp);
			
			printf("done\n");
		}
	
		void save_track_float(string filename){
				
			printf("save_track_float('%s')...", filename.c_str()); fflush(stdout);
			
			FILE *fp = _fopen(filename, "wb");
			fwrite(track, sizeof(*track), n/tracks_freq*ns/track_freq*np*4, fp);	
			fclose(fp);
			
			printf("done\n");
		}
	
		void save_gadget1(string filename) {
				
			GF1IO gsio;
				
			gsio.header.time = 0;
			gsio.header.npart[0] = gsio.header.npartTotal[0] = 0;
			gsio.header.npart[1] = gsio.header.npartTotal[1] = halo.n;
			gsio.header.npart[2] = gsio.header.npartTotal[2] = disk.n;
			gsio.header.npart[3] = gsio.header.npartTotal[3] = 0;
			gsio.header.npart[4] = gsio.header.npartTotal[4] = 0;
			gsio.header.npart[5] = gsio.header.npartTotal[5] = 0;
	
			gsio.header.mass[0] = 0;
			gsio.header.mass[1] = halo.n ? halo.m/halo.n : 0; 
			gsio.header.mass[2] = disk.n ? disk.m/disk.n : 0;
			gsio.header.mass[3] = 0;
			gsio.header.mass[4] = 0;
			gsio.header.mass[5] = 0;

			gsio.header.npartTotalHighWord[0] = 0;
			gsio.header.npartTotalHighWord[1] = 0;
			gsio.header.npartTotalHighWord[2] = 0;
			gsio.header.npartTotalHighWord[3] = 0;
			gsio.header.npartTotalHighWord[4] = 0;
			gsio.header.npartTotalHighWord[5] = 0;
				
						
			gsio.header.flag_stellarage = 0;	
			gsio.header.flag_metals = 0;
		
			gsio.header.flag_entropy_instead_u = 0;	
			gsio.header.flag_doubleprecision = 0;	

			gsio.header.flag_ic_info = 0;          
			gsio.header.lpt_scalingfactor = 0;   
				
			gsio.header.num_files=1;
				
			gsio.NumPart = halo.n+disk.n;
			gsio.Ngas = 0;
			gsio.N = gsio.NumPart;

			gsio.print();

			gsio.p[0] = new GF1IO::P[gsio.NumPart];
			gsio.Id = new uint[gsio.NumPart];

			if (halo.n) gsio.p[1] = &gsio.p[0][0];
			if (disk.n) gsio.p[2] = &gsio.p[0][halo.n];

			// halo
			for (uint i=0; i<halo.n; i++) {
					
				for (uint k=0; k<3; k++) {
					gsio.p[1][i].Pos[k] = halo.p[i].Pos[k];
					gsio.p[1][i].Vel[k] = halo.p[i].Vel[k]==halo.p[i].Vel[k] ? halo.p[i].Vel[k] : 0;
				}
	
				gsio.p[1][i].Type = 1;
				gsio.Id[i] = i+1;
			}
				
			// disk
			for (uint i=0; i<disk.n; i++) {
				for (uint k=0; k<3; k++) {
					gsio.p[2][i].Pos[k] = disk.p[i].Pos[k];
					gsio.p[2][i].Vel[k] = disk.p[i].Vel[k]==disk.p[i].Vel[k] ? disk.p[i].Vel[k] : 0;
				}
				gsio.p[2][i].Type = 2;
				gsio.Id[halo.n+i] = halo.n+i+1;
			}
				
				
			gsio.print_P(15);
			gsio.write(filename.c_str());
		
		}
	
		void load_gadget1(const char *fname) {
			
			printf("load_gadget1('%s')...", fname); fflush(stdout);
			
			
			GF1IO gsio;
			gsio.read(fname);
			
			halo.free();
			halo.n = gsio.header.npart[1];
			halo.p = new P[halo.n];
			halo.m = gsio.header.mass[1]*halo.n;
			
			for (uint i=0; i<halo.n; i++)
				for (uint k=0; k<3; k++) {
					halo.p[i].Pos[k] = gsio.p[1][i].Pos[k];
					halo.p[i].Vel[k] = gsio.p[1][i].Vel[k];
					halo.p[i].Mass = gsio.header.mass[1];
				}
			
			disk.free();
			disk.n = gsio.header.npart[2];
			disk.p = new P[disk.n];
			disk.m = gsio.header.mass[2]*disk.n;
			
			for (uint i=0; i<disk.n; i++)
				for (uint k=0; k<3; k++) {
					disk.p[i].Pos[k] = gsio.p[2][i].Pos[k];
					disk.p[i].Vel[k] = gsio.p[2][i].Vel[k];
					disk.p[i].Mass = gsio.header.mass[2];
				}
			
			printf("GF1IO::Mtot =%g\n", gsio.header.mass[1]*gsio.header.npart[1] + gsio.header.mass[2]*gsio.header.npart[2]);
		}
	
		void saveAMRGrid(const char* fname, vector< vector<int> > &map, REAL *m){
			
			float *_m = new float[dg.n];
			
			for (int i=0; i<map.size(); i++ ) {
				double mi(0);
				for (int j=0; j<map[i].size(); j++ )  mi+= m[map[i][j]];
				for (int j=0; j<map[i].size(); j++ ) _m[map[i][j]] = mi;
			}
			
			save_arr(fname, _m, dg.n);
			delete [] _m;
		}
		
		
		
		
		void print() {
			printf("\n");
			printf("ICs::print()\n");
			printf(" ├ n = %i\n", n);
			printf(" ├ ns = %i\n", ns);
			printf(" ├ np = %i\n", np);
			printf(" │\n");
			printf(" ├ dg.nr = %i\n", dg.nr);
			printf(" ├ dg.nz = %i\n", dg.nz);
			printf(" │\n");
			printf(" ├ dvk_mu = %g\n", dvk_mu);
			printf(" ├ dvh_mu = %g\n", dvh_mu);
			printf(" │\n");
			printf(" ├ dh_min = %g\n", dh_min);
			printf(" ├ dh_max = %g\n", dh_max);
			printf(" │\n");
			printf(" ├ step_out = %i\n", step_out);
			printf(" ├ h = %g\n", h);
			printf(" ├ stype = %i\n", stype);			
			printf(" ├ curfolder = %s\n", curfolder);
			printf(" ├ ssfolder = %s\n", ssfolder);
			printf(" └ nsteps = %i\n", nsteps);
			printf("\n");
		}
	

		void save(int i) {
			
			char fname[256];
			sprintf(fname, "dump.%05i.dat", i);
			printf("save(%s)...", fname);  fflush(stdout);
			
			FILE *fp = _fopen(fname, "wb");
			fwrite(this, sizeof(*this), 1, fp);
			fwrite(p, sizeof(*p), n, fp);
			
			fclose(fp);
			
			if (i) {
				sprintf(fname, "%s/dump.%05i.dat", curfolder, i-1);
				remove(fname);
			}
			
			printf("done\n");
		};
		
		void load(const char *fname) {
			
			printf("load(%s)...", fname);  fflush(stdout);
			
			FILE *fp = _fopen(fname, "rb");
			fread(this, sizeof(*this), 1, fp);
			
			p = new P[n];
			fread(p, sizeof(*p), n, fp);
			
			fclose(fp);
			
			printf("done\n");
			
		};
	
	private:

		FILE* _fopen(string filename, const char* mode){
			filename = string(curfolder) + "/" + filename;
			string path = filename.substr(0, filename.find_last_of("/"));
			mkdir(path.c_str(), 0777);
			return fopen(filename.c_str(), mode);
		}
		
		
};


#endif
